function L = LEs(alpha, beta, xi, mu, A, omega, dt, T0, Tsim)
% LEs  计算最大李雅普诺夫指数
%   L = LEs(alpha, beta, xi, mu, A, omega, dt, T0, Tsim)
%   输入：
%     alpha, beta, xi, mu, A, omega — 模型参数
%     dt    — 时间步长
%     T0    — 丢弃瞬态时间
%     Tsim  — 统计时间
%   输出：
%     L — 最大李雅普诺夫指数

% 计算步数
N0    = floor(T0/dt);
Nsim  = floor(Tsim/dt);

% 初始状态与扰动向量
S = [0.2; 0.1; 0.2];
v = randn(3,1); v = v/norm(v);

sumLog = 0;
count  = 0;

% 每步正交化
for k = 1 : (N0 + Nsim)
    t = (k-1)*dt;
    mu_s = A * cos(omega*t);
    
    % —— 主系统和变分——RK4耦合积分 —— 
    % k1
    k1S = f(S, alpha, beta, xi, mu, mu_s);
    J1  = Jacobian(S, alpha, beta, xi, mu, mu_s);
    k1v = J1 * v;
    % k2
    S2  = S + 0.5*dt*k1S;
    v2  = v + 0.5*dt*k1v;
    t2  = t + 0.5*dt; mu_s2 = A*cos(omega*t2);
    k2S = f(S2, alpha, beta, xi, mu, mu_s2);
    J2  = Jacobian(S2, alpha, beta, xi, mu, mu_s2);
    k2v = J2 * v2;
    % k3
    S3  = S + 0.5*dt*k2S;
    v3  = v + 0.5*dt*k2v;
    t3  = t + 0.5*dt; mu_s3 = A*cos(omega*t3);
    k3S = f(S3, alpha, beta, xi, mu, mu_s3);
    J3  = Jacobian(S3, alpha, beta, xi, mu, mu_s3);
    k3v = J3 * v3;
    % k4
    S4  = S +     dt*k3S;
    v4  = v +     dt*k3v;
    t4  = t +     dt; mu_s4 = A*cos(omega*t4);
    k4S = f(S4, alpha, beta, xi, mu, mu_s4);
    J4  = Jacobian(S4, alpha, beta, xi, mu, mu_s4);
    k4v = J4 * v4;
    
    % 更新
    S = S + (dt/6)*(k1S + 2*k2S + 2*k3S + k4S);
    v = v + (dt/6)*(k1v + 2*k2v + 2*k3v + k4v);
    
    % 丢弃瞬态后每步正交化并累积
    if k > N0
        nv      = norm(v);
        sumLog  = sumLog + log(nv);
        v       = v / nv;
        count   = count + 1;
    end
end

% 输出最大李雅普诺夫指数
L = sumLog / (count * dt);

end

